Monte  Carlo  Modeling  of  YBCO  Vapor  Deposition 


Jing  Fan  ,  Iain  D.  Boyd",  and  Chris  Shelton' 


Department  of  Aerospace  Engineering,  University  of  Michigan,  Ann  Arbor,  MI  48109 
3M  Corporate  Process  Technology,  St.  Paul,  MN  551 44* 


Abstract.  A  three  dimensional  electron  beam  physical  vapor  deposition  process  of  super-conducting  films  of 
YBayCu^O-y-g  in  a  vacuum  chamber  is  investigated  both  computationally  and  experimentally.  The  numerical  analysis 
employs  the  direct  simulation  Monte  Carlo  (DSMC)  method.  Some  important  modeling  issues  such  as  atomic  collision 
cross-sections  for  metal  vapors  and  hyperfine  electronic  structure  of  the  atomic  absorption  spectra  are  addressed.  Film 
deposition  thickness  profiles  and  atomic  absorption  spectra  given  by  the  DSMC  method  and  experiment  for  pure  yttrium 
evaporation  are  in  excellent  agreement.  Interaction  between  the  vapor  jets  is  found  to  have  a  significant  effect  on  the  film 
growth  rate  and  species  molar  ratio  that  are  key  factors  of  the  film  technology. 


INTRODUCTION 


Electron  Beam  Physical  Vapor  Deposition  (EBPVD)  is  an  important  deposition  technique  for  a  variety  of 
materials  that  involve  the  use  of  an  electron  beam  to  vaporize  metallic  atoms  from  a  solid  ingot.  Applications  of  this 
technology  include  coating  of  high  temperature  super-conducting  films,  metallic  coating  of  fibers,  and  enrichment  of 
uranium  for  nuclear  power  generation.  A  critical  component  of  these  processes  is  the  manner  in  which  the  atoms  are 
transported  from  the  ingot  to  the  deposition  surface.  There  is  a  need  to  understand  in  detail  the  gas  dynamics  of  the 
expansion  process  particularly  in  cases  where  multiple  ingots  of  different  species  are  evaporated. 

In  this  paper,  we  consider  both  numerically  and  experimentally  the  deposition  of  super-conducting  films  of 
YBa2CujOy_g  that  are  of  great  technological  interest  in  numerous  applications.1  Figure  1  shows  the  top  and  side 
views  of  the  experimental  deposition  chamber.  Using  vacuum  pumps,  the  background  pressure  is  maintained  at 

2  ~  5xlO~sTorr .  In  the  EBPVD  process,  atoms  are  vaporized  from  the  molten  surfaces  of  yttrium,  barium  and 
copper  ingots  by  a  high-energy  electron  beam.  The  vapor  jets  expand  through  the  chamber,  impinge  and  deposit  on  a 
substrate.  The  vaporization  rates  achieved  with  the  electron  beam  produce  flow  at  relatively  high  density  close  to  the 
source  surfaces.  However,  as  the  jets  proceed  toward  the  substrate,  rapid  expansion  brings  the  flow  into  the  non¬ 
continuum  regime.  The  direct  simulation  Monte  Carlo  (DSMC)  method2  is  used  to  simulate  the  three-dimensional 
process.  Besides  the  translational  energy  mode,  the  atomic  electronic  energy  is  taken  into  account.3  At  this  stage,  no 
flux  of  oxygen  is  included. 
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FIGURE  1 .  Schematic  diagrams  of  the  deposition  chamber.  Left :  top  view;  right:  side  view. 
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ATOMIC  COLLISION  CROSS-SECTIONS  FOR  METAL  VAPORS 


The  DSMC  method  is  a  powerful  and  general  numerical  technique  for  computation  of  low  density, 
nonequilibrium  flows,  and  has  been  developed  over  the  past  30  years  primarily  for  aerospace  applications.2 
Application  of  the  DSMC  method  to  metal  vapor  flows  is  straightforward  in  principle  and  has  been  demonstrated  in 
relatively  simple  geometry  such  as  an  axisymmetric  case.'1  Some  specific  complications  do  naturally  arise.  One  of 
them  is  the  determination  of  atomic  collision  cross-sections  for  metal  vapors.  Appropriate  collision  cross-sections 
are  essential  for  the  DSMC  method  to  simulate  particle  collisions.  They  are  usually  derived  from  coefficients  of 
viscosity.2,4  The  measurement  of  the  transport  properties  for  metal  vapors  at  high  temperature,  however,  has  proven 
very  difficult  because  of  many  unavoidable  sources  of  error.5  Actually,  except  for  alkali  metal  vapors,610  no 
viscosity  data  are  available  for  metal  vapors  such  as  yttrium,  barium  and  copper,  and  few  studies  theoretically 
address  the  collision  cross-sections  for  metal  vapors. 

The  present  paper  employs  the  inverse-power  model  to  describe  the  interaction  between  atoms  of  metal  vapors 


F  =  k/^  . 

The  collision  cross-section  and  coefficient  of  viscosity  for  the  inverse-power  model  may  be  written  as4 
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where  W0m  is  a  constant,  mr  is  the  reduced  mass,  cr  is  the  relative  speed,  k  is  the  Boltzmann  constant,  mis  the 


atomic  mass,  T  is  the  temperature,  and  A2(tj )  is  a  numerical  factor. 

The  values  of  rj  and  K  for  alkali  metal  vapors  have  been  obtained  by  equating  the  right-hand  side  of  Eq.  (3)  with 
theoretical  data10  at  temperatures  of  700K  and  2000K,  respectively,  and  they  are  given  in  Table  1.  Figure  2  compares 
the  profiles  of  Eq.  (3)  to  the  theoretical  and  experimental  data  for  sodium  and  cesium.  The  experimental  data,3 
though  with  scatter,  support  in  trend  the  theoretical  data10  and  fitted  profiles  that  are  in  excellent  agreement. 


TABLE  1.  Parameters  of  the  inverse-power  model  for  alkali  metal  vapors. 


Element 

n 

K 

A2(il) 

Lithium 

6.333 

7.984xlO~70 

0.379 

Sodium 

7.351 

4.979xl0~79 

0.355 

Potassium 

6.737 

4.462xl0~73 

0.368 

Rubidium 

6.669 

2.521xl0~72 

0.370 

Cesium 

6.531 

7. 081  x  IQ-71 

0.374 

FIGURE  2.  Comparison  of  theoretical,  calculated  and  experimental  viscosity  for  sodium  and  cesium  vapors. 
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DSMC  calculations  often  employ  the  variable  hard-sphere  (VHS)  model2 11  to  describe  the  interaction  between 
particles.  It  assumes  that  the  scattering  during  particle  collisions  is  isotropic  in  the  center  of  mass  frame  of  reference, 
but  the  collision  cross-section  has  the  same  relation  with  the  relative  speed  as  the  inverse-power  model 
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where  ,uref  is  the  coefficient  of  viscosity  at  a  reference  temperature Tref  ,  and  w  =  0.5  +  2/(jj-1  ) . 

The  reference  collision  diameters  drej ■  for  alkali  metal  vapors  have  been  obtained  with  /Jrej  calculated  using  Eq. 
(3),  and  they  are  given  in  Table  2.  We  assume  that  the  values  of  //and  K for  the  alkali  metal  vapors  apply  to  other 
metal  vapors  occupying  the  same  row  in  the  periodic  table  of  elements.  Based  on  the  assumption,  the  reference 
collision  diameters  for  several  metal  vapors  of  interest  are  calculated  and  given  in  Table  2.  The  selection  of  Trer  is 

arbitrary  and  does  not  affect  collision  cross-sections.  To  be  reasonable  physically,  it  is  taken  to  be  2,0()()K  that  is 
above  the  melting  point  of  all  elements  considered  in  Table  2. 

An  outstanding  feature  of  the  reference  collision  diameters  for  the  metal  vapors  is  that  they  are  quite  large  in 
comparison  with  common  atomic  gases.  This  is  mainly  due  to  the  larger  viscosity  of  the  gases.  For  example,  the 

coefficient  of  viscosity  for  argon  gas  is  5.452x10  '  Nms  at  700K  and  O.IMPa,  while  the  viscosity  of  sodium 
vapor  at  the  same  conditions  is  1.475xl0~5  Nms~2  ,10 


TABLE  2.  Parameters  of  the  VHS  model  for  metal  vapors  at  Tref  =  2000 K  . 


Element 

dref  xlO  10  m 

(0 

Lithium 

4.562 

0.875 

Sodium 

4.961 

0.815 

Potassium 

5.554 

0.849 

Rubidium 

5.846 

0.853 

Cesium 

6.306 

0.862 

Copper 

6.272 

0.849 

Titanium 

5.844 

0.849 

Yttrium 

5.903 

0.853 

Zirconium 

5.941 

0.853 

Barium 

6.358 

0.862 

ATOMIC  ABSORPTION  SPECTRA 

In  situ  monitoring  and  control  is  one  of  the  primary  issues  involved  in  EBPVD  processing  of  super-conducting 
films  of  YB(12Cut)0'i^§  .  Diode-laser-based  absorption  monitors  for  yttrium  and  barium  have  been  developed  to 
meet  this  need.  13  They  are  simultaneously  both  element-specific  and  non-invasive,  and  can  provide  the  atomic 
absorption  spectra  from  which  the  atomic  flux  across  the  laser  path  may  be  derived. 

Using  the  monitors,  an  atomic  absorption  spectral  distribution  around  a  central  frequency  v0  along  a  laser  path  is 

observed.  The  spectrum  I(v )  in  the  vicinity  of  v(j  is  called  the  line  profile.14  The  profile  may  have  many  peaks 
referred  to  as  hyperfine  electronic  structure.  They  correspond  to  the  frequencies  absorbed  by  particles  when  the  light 
beam  traverses  the  subject  material,  and  are  centered  at  Vj  =  Vq  +  Avt  ( Avt  «  v() ).  These  peaks  have  a  finite 
•  _2 

width  due  to  various  effects  of  broadening.  In  low  pressure  (10  Torr )  situations,  the  width  is  mainly  determined 
by  the  Doppler  effect.14  This  takes  place  because  frequencies  absorbed  by  particles  with  different  velocities  in  the 
light  beam  direction  are  slightly  different. 

At  thermal  equilibrium,  a  Doppler-broadening  peak  shape  centered  at  a  frequency  vt  is  described  by14 
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(5) 


Ii(v)  =  ^exp[-p2Z2i(v- vf  )2  ] , 
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where  the  wavelength  X-t  =  c/vt  ,  c  is  the  light  speed,  and  /3  =  -Jm/YkT  ,  The  spectral  distribution  (5)  corresponds 
to  a  Maxwellian  velocity  distribution 


f(ve  )  =  -j=  exp{p2  v2e), 

Vtt 


(6) 


where  is  the  thermal  velocity  component  in  the  light-beam  direction.  Comparing  the  distributions  (5)  and  (6) 
leads  to  the  relations 


V-Vi  =V(J?ii  , 


(7) 


and 


Ii(v)  =  Aif(ve).  (8) 

A  velocity  distribution  of  particles  may  be  obtained  from  DSMC  calculations.  To  compare  with  an  experimental 
absorption  spectrum,  it  needs  to  be  converted  into  an  absorption  scan  using  equations  (7)  and  (8).  Before  the 
conversion,  we  must  normalize  the  calculated  velocity  distribution  f(  vg  )  in  accordance  with  the  measured 
absorption  spectrum  /,(  v ) .  The  normalization  condition  requires 

An  bVf>  Vf-  )v*f-  'dve=1’  (9) 

and 

Bv\avl(v)v~‘ dv  =  l ,  (10) 

where  vj  and  v  are  velocity  and  frequency  normalization  factors,  respectively,  and  Avf  and  Bv  are  overall 
normalization  factors.  The  relation  (7)  requires 

v'=v;A.  (id 

To  compare  calculated  and  experimental  atomic  absorption  spectra  properly,  hyperfine  electronic  structure  must 
be  included.14,  15  For  example,  an  experimental  profile  line  for  the  ground  state  transition  of  yttrium  at 

Vy  =4. 486x10 14  Hz  consists  of  four  peak  components  with  Av-t  and  relative  intensity  factors  Yi  given  in  the 
second  column  of  Table  3.  Accordingly,  the  calculated  line  profile  must  also  be  the  weighted  sum  of  the  lour  peak 
shapes  centered  at  slightly  different  frequencies  v,  =  Vy  +  AVj 
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cal(v)=  I Yilfa'(v ) 
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TABLE  3.  Hyperfine  structure  parameters  of  atomic  absorption  spectra  for  yttrium. 


VyXlO14  (Hz) 

4.486 

4.416 

Av,xl08  (Hz) 

0.00 

0.00 

Av2xl08  (Hz) 

1.14 

0.86 

Av3xlOS  (Hz) 

1.76 

1.99 

Av4x108  (Hz) 

2.90 

1.99 

Yi 

5 

20 

Yi 

1 

1 

Yi 

1 

1 

Y4 

9 

4 
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RESULTS  AND  DISCUSSION 


Two  cases  are  considered  here.  For  the  first  case  where  only  the  yttrium  source  is  evaporating,  the  evaporation 
rates  are  measured  by  the  feed  rate  of  the  ingot  into  the  chamber  (Table  4).  The  source  surface  temperatures  Ts  are 
calculated  from16 

log I()  rg=C -0.5 log  10  Ts  - B/Ts  ,  (13) 

where  r„  denotes  the  evaporation  rate  in  unit  of  gj cm2  sec .  The  values  of  B  and  C  for  many  metal  vapors 
including  yttrium,  barium  and  copper  are  given  in  Table  10.2  in  Ref.  16.  The  number  densities  at  the  source 
surfaces  are  calculated  from  n  =  (3re ,  where  fe  =  rg  /m ,  and  ji  =  y Jm  /  2kTs  . 


2 

TABLE  4.  Source  evaporation  rates  re(  moles  /  cm  sec)  and  surface  temperature  Ts  (  K  ) . 


Case 

Yttrium 

Barium 

Copper 

re 

T 

1  s 

r  T 

1  e  1  s 

re 

T 

1  s 

I 

1.129xl0~5 

2091 

0  — 

0 

— 

II 

8.542xl0~6 

2065 

1.708xl0~5  1010 

2.562xl0~5 

1732 

DSMC  calculations  are  carried  out  using  a  software  system  called  MONACO- YBC. 17  In  an  orthogonal 
coordinate  system  with  the  origin  at  O,  the  x-  and  y-  axes  are  along  OO'  and  OA  respectively,  and  the  z  axis  points 
to  the  substrate  (cf.  Fig.  1).  The  computational  grid  employs  20x60x50  cuboid  cells.  The  cell  sizes  are  non- 
uniform:  cell  size  in  the  z  direction  near  the  source  surfaces  are  very  small  to  capture  the  rapid  jet  expansion,  and 
cell  sizes  in  the  x  and  y  directions  in  the  region  covering  the  sources  are  small  to  resolve  the  region  of  flow  that  is 
most  important  for  accurate  simulation  of  the  deposition  process.  About  two  million  simulated  particles  are  used  in 
the  calculations.  The  results  may  be  obtained  overnight  on  a  SGI  OCTANE  workstation  or  a  DELL/Pentium  III  PC. 

The  chamber  walls  and  block  plate  surfaces  are  assumed  to  be  perfectly  sticking.  This  assumption  is  well 
satisfied  in  this  study  because  the  wall  and  plate  surface  temperatures  are  low.  New  particles  are  vaporized  into  the 
chamber  from  the  surfaces  of  the  sources.  They  may  collide  with  one  another  above  the  pool,  particularly  for  cases 
with  high  evaporation  rates.  This  causes  a  fraction  of  the  atoms  to  develop  post-collision  velocities  that  move  them 
back  to  the  molten  pools.  In  Ref.  3,  a  total  back  scatter  flux  was  found  to  be  as  high  as  19%  of  the  evaporation  rate. 
To  maintain  a  net  inflow  rate  measured  by  the  feed  rate  of  an  ingot  into  the  chamber,  the  back-scattered  particles  are 
made  to  diffusely  reflect  from  the  source  surface. 

Case  I:  Evaporation  of  Pure  Yttrium 

An  experiment  of  pure  yttrium  with  evaporation  rate  given  in  Table  4  is  carried  out  that  lasts  12  minutes.  Figure 
3  compares  measured  deposition  thickness  data  with  DSMC  and  collisionless  profiles.  The  DSMC  and  collisionless 
calculations  are  carried  out  using  the  same  code,  but  the  latter  turns  off  atomic  collisions.  The  deposition  thickness  is 
measured  with  a  profilometer  along  a  symmetrical  line  of  the  substrate.  The  substrate  used  is  a  polished  silicon 
wafer.  Since  the  coating  is  rather  thick  and  blistering  of  the  film  does  not  allow  accurate  sampling  points,  scatter  in 
the  measured  data  is  clearly  seen.  Deducting  the  measured  data  scatter,  the  DSMC  profile  compares  well  with  the 
experiment.  The  collisionless  profile  is  significantly  lower  than  the  DSMC  and  measured  results.  The  atomic 
collisions  impede  to  some  extent  the  atomic  diffusion  to  the  block  plates  and  chamber  walls  that  are  assumed  to  be 
perfectly  sticking.  Neglecting  these  collisions  results  in  the  atoms  having  more  chances  to  stick  on  the  block  plate 
surfaces  and  chamber  walls  than  to  be  deposited  on  the  substrate. 

Figure  4  shows  DSMC  and  measured  atomic  absorption  (AA)  spectra  across  an  aperture  close  to  the  symmetric 
line  of  the  substrate.  Because  the  hyperfine  electronic  structures  of  yttrium  are  different  at  the  frequencies  of  668nm 
and  679nm  (Table  3),  the  corresponding  AA  peak  structures  are  inclined  to  the  right  and  left,  respectively.  The 
DSMC  and  measured  Doppler  widths  and  peak-structure  details  are  all  in  excellent  agreement. 
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Figure  5  shows  the  number  density  given  by  the  DSMC  calculation.  It  shows  that  there  is  a  rapid  expansion 
around  the  source  that  leads  to  a  significant  decrease  of  the  number  density.  Two  continuum  breakdown  criteria, 
proposed  by  Bird18  and  Boyd,  et  al.19,  respectively,  are  computed  based  on  the  DSMC  number  density  field.  For 
steady  flows  Bird’s  criterion  may  be  written  as 


P  =  Ma. 


(14) 


where  Ma,X,p  are  the  local  Mach  number,  mean  free  path,  and  density,  respectively,  y  is  the  specific  heat  ratio. 


and  s  is  the  distance  along  a  streamline.  Studies  in  Ref.  18  indicated  that  the  continuum  approach  in  expanding  flows 
failed  at  a  value  for  P  of  about  0.05.  Comparing  Navier-Stokes  and  DSMC  results  for  shock  waves  and  hypersonic 
flows  around  a  sphere,  Boyd,  Chen  and  Candler19  found  that  the  gradient-length  local  Knudsen  number  based  on 
density 

AlVpl 

^nGLL—D  =  —  ’  (15) 

P 


can  more  precisely  specify  the  limit  at  which  the  Navier-Stokes  equations  were  invalid.  The  continuum  approach 
broke  down  wherever  the  value  of  KnGLL_D  exceeded  0.05.  19  Variation  of  the  two  criteria  along  a  path  from  the 
yttrium  source  center  perpendicularly  to  the  substrate  is  showed  in  Fig.  5.  Both  of  them  indicate  that  the  flow  is 
completely  in  the  non-continuum  regime  that  needs  to  be  studied  kinetically. 


FIGURE  3.  Film  deposition  thickness  distributions  for  Case  I.  Left:  comparison  of  DSMC,  collisionless  and  measured  results 
along  a  symmetrical  line  of  the  substrate;  right:  distribution  ( flm  )  over  the  whole  substrate  given  by  the  DSMC  calculation. 


v-Vj  (Hz)  v-v„  (Hz) 

FIGURE  4.  Comparison  of  DSMC  (solid  line)  and  measured  (triangle)  atomic  absorption  spectra  for  Case  I.  Left:  at  frequency 
of  4.486 XlO 14  Hz  ( 668 nm  ) ;  right:  at  frequency  of  4.416  XlO 14  Hz  ( 679nm  )  . 
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FIGURE  5.  Left :  DSMC  number  density  field  (m'3)  for  Case  I;  right :  variation  of  two  continuum  breakdown  criteria,  derived 
from  the  left  number  density  field,  along  a  path  from  the  yttrium  source  center  perpendicularly  to  the  substrate. 


Case  II:  Co-evaporation  of  Yttrium,  Barium  and  Copper 

The  number  densities  of  the  yttrium,  barium  and  cooper  sources  computed  from  the  evaporation  rates  given  in 

Table  4  are  2.903xl020 ,  1.032xl021  and  8.401xl020 m~3 ,  respectively.  Figure  6  (left)  shows  the  total  number 
density  field  contours.  Similar  to  Case  I,  there  is  a  significant  decrease  of  the  number  density  around  the  sources  that 
corresponds  to  a  rapid  expansion. 


FIGURE  6.  DSMC  results  for  Case  II.  Left:  total  number  density  field  (m 3);  right:  comparison  of  species  molar 
ratio  along  a  symmetrical  line  of  the  substrate.  Collision  — ,  part-collision  — ,  and  collisionless  . 


To  understand  the  effect  of  collisions  between  different  species,  molar  species  ratios  of  yttrium,  barium,  and 
copper  across  a  symmetrical  line  of  the  substrate  are  compared  at  collision,  part-collision  (turning  off  collisions 
between  the  different  species),  and  collisionless  situations  (Fig.  6,  right).  The  part-collision  profiles  are  close  to  the 
collisionless  profiles  that  are  relatively  uniform  and  keep  a  ratio  near  to  the  original  one  between  the  yttrium  barium, 
and  copper  sources,  i.e.,  1:  2:  3.  Because  the  collisions  between  yttrium  and  barium  or  copper  atoms  generally 
impede  a  diffusion  of  the  yttrium  atoms  in  the  positive  y-direction,  the  collision  profile  of  yttrium  is  significantly 
non-uniform  and  peaks  at  y  =  0  .  A  similar  mechanism  results  in  a  peak  at  y  =  18”  in  the  collision  profile  of  copper. 
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Since  the  number  density  of  the  copper  source  is  about  three  times  that  of  the  yttrium,  the  collisions  between  barium 
and  copper  atoms  are  more  frequent  than  those  between  barium  and  yttrium  atoms.  This  causes  the  barium  peak  to 
be  inclined  to  the  left.  The  significantly  non-uniform  distributions  of  the  molar  ratio  are  actually  deleterious  to  the 
quality  of  super-conducting  films  of  YBa2Ciij07_g that  requires  a  ratio  of  1:2:3  between  yttrium,  barium  and 
copper.  The  interaction  between  the  vapor  jets  revealed  by  the  calculations  is  helpful  to  improve  the  film  technology. 


CONCLUSIONS 

A  three  dimensional  electron  beam  physical  deposition  of  super-conducting  films  of  YBa2Cu707_g  was  studied 
using  the  DSMC  method  and  experiment.  Atomic  collision  cross-sections  for  metal  vapors  such  as  yttrium,  barium 
and  copper  were  given.  They  were  quite  large  in  comparison  with  those  for  common  atomic  gases.  A  means  to 
convert  a  velocity  distribution  of  atoms  into  an  absorption  spectrum  was  established.  Film  deposition  thickness 
profiles  along  a  symmetrical  line  of  the  substrate  and  atomic  absorption  spectra  across  an  aperture  close  to  the 
symmetrical  line  given  by  the  DSMC  method  and  experiment  for  pure  yttrium  evaporation  were  in  excellent 
agreement.  Interaction  between  the  vapor  jets  was  found  to  have  a  significant  effect  on  the  film  growth  rate  and 
species  molar  ratio  that  are  key  factors  of  the  film  technology. 
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